#
# example028.py
#
# Copyright (C) 2012 Robert Buj Gelonch
# Copyright (C) 2012 David Megias Jimenez
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
__author__ = "Robert Buj Gelonch, and David Megias Jimenez"
__copyright__ = "Copyright 2012, Robert Buj Gelonch and David Megias Jimenez"
__credits__ = ["Robert Buj Gelonch", "David Megias Jimenez"]
__license__ = "GPL"
__version__ = "3"
__maintainer__ = "Robert Buj"
__email__ = "rbuj@uoc.edu"
__status__ = "Development"
__docformat__ = 'plaintext'

from numpy import arange
from numpy import array
from numpy import concatenate
from numpy import cos
from numpy import dot
from numpy import pi
from numpy import sin
from numpy import zeros
from pylab import figure
from pylab import grid
from pylab import plot
from pylab import show
from pylab import title
from pylab import xlabel
from scipy.linalg import toeplitz
from scipy.signal import firwin
from scipy.signal import kaiserord

#----------------------------------------------------------
# Return the filtered signal convolving the impulsive res-
# ponse with the input signal by using the toeplitz matrix
#----------------------------------------------------------
def concatenate_filtering(taps, x):
    m = len(taps) + len(x)-1
    hc = concatenate((taps, zeros(len(x))))
    hr = concatenate(([taps[0]], zeros(m)))
    h = toeplitz(hc, r=hr)
    y = dot(h, array([concatenate((x, zeros(len(taps))))]).transpose())
    return y

#------------------------------------------------
# Create a signal for demonstration.
#------------------------------------------------
sample_rate = 100.0
nsamples = 400
t = arange(nsamples) / sample_rate
x = cos(2 * pi * 0.5 * t) + 0.2 * sin(2 * pi * 2.5 * t + 0.1) + \
    0.2 * sin(2 * pi * 15.3 * t) + 0.1 * sin(2 * pi * 16.7 * t + 0.1) + \
        0.1 * sin(2 * pi * 23.45 * t + .8)

#------------------------------------------------
# Create a FIR filter and apply it to x.
#------------------------------------------------
# The Nyquist rate of the signal.
nyq_rate = sample_rate / 2.0
# The desired width of the transition from pass to stop,
# relative to the Nyquist rate.  We'll design the filter
# with a 5 Hz transition width.
width = 5.0 / nyq_rate
# The desired attenuation in the stop band, in dB.
ripple_db = 60.0
# Compute the order and Kaiser parameter for the FIR filter.
N, beta = kaiserord(ripple_db, width)
# The cutoff frequency of the filter.
cutoff_hz = 10.0
# Use firwin with a Kaiser window to create a lowpass FIR filter.
taps = firwin(N, cutoff_hz / nyq_rate, window=('kaiser', beta))
# Apply the FIR filter to the input signal.
y = concatenate_filtering(taps, x)

#------------------------------------------------
# Plot the original and filtered signals.
#------------------------------------------------
# The phase delay of the filtered signal.
delay = 0.5 * (N-1) / sample_rate
f1 = figure(1)
# Plot the original signal.
plot(t, x)
# Plot the filtered signal, shifted to compensate for the phase delay.
t_y = arange(len(y)) / sample_rate
plot(t_y-delay, y, 'r-')
# Plot just the "good" part of the filtered signal.  The first N-1
# samples are "corrupted" by the initial conditions.
plot(t_y[N-1:-(N-1)]-delay, y[N-1:-(N-1)], 'g', linewidth=4)

title('Digital filtering in time domain: toeplitz matrix')
xlabel('t')
grid(True)
f1.savefig("../plots/example028.eps", format="eps")

show()